The Discrete Element Method (DEM) is a process executed by a computer system that allows the prediction of the flow and structure of granular (particulate) materials. The principle physics associated with the process is that of collisions between particles and between particles and solid bodies (such as equipment). Examples of applications of the process include, but are not limited to:                (i) Crushing and grinding of rocks in mineral processing, where the data produced by the process provides information about the flow and breakage of particles for the purposes of improving process efficiency, production rate and the equipment design and reducing equipment wear;        (ii) Mixing of powders and grains for food, plastics, pharmaceuticals and household product manufacture, where the data produced by the process provides information about flow, mixing, particle growth, breakage and agglomeration for the purposes of improving the mixing process and the equipment design;        (iii) Transport and flow of particles through processing equipment, where the data produced by the process provides information used to understand and optimise equipment design;        (iv) Digging and excavation, where the data produced by the process provides information used to improve efficiency;        (v) Rock cutting and fragmentation, where the information is used to improve cutting efficiency and reduce equipment wear;        (vi) Landslides and debris flow, where the data produced by the process provides information used for planning, mitigation and disaster response management.        
In a DEM simulation there will be many particles, each with specific attributes such as size, shape, density and collision properties. There will also often be boundary objects which can have an arbitrarily complex shape typically represented as surface meshes or by trimmed geometric primitives such as lines, planes, cylinders, arcs etc.
A basic DEM process proceeds by time marching, where the system evolves from one known state to new one by a small increment in time, and an update of all the state information is generated. There are three elements:                (i) determine which particles are in contact;        (ii) for particles that are in contact with other particles or boundaries calculate the forces between them; and        (iii) sum the net forces and torques on the particles and objects and integrate their motion to predict new positions, orientations, velocities and spins.        
The process for determining which particles are in contact can be an enormously complex and computationally expensive part of the method and so specific approaches have been developed to reduce these costs to manageable levels.
The simplest approach, whose cost is of the order of N2, where N is the number of particles, is to simply check each pair of possible contact pairs at each timestep. For small values of N (say <1000) this is feasible but very inefficient.
For larger simulations the cost of the contact detection needs to be of order N for the simulation to be feasible. This may be achieved by performing the search for particles that are in contact less frequently and during the search by building a list of pairs of entities (particles or parts of boundary objects) that may collide during a particular period of time from the time of the search. This is sometimes termed the neighbour list or near neighbour list or interaction list. At each timestep the contact detection is performed only for pairs of entities in the neighbour list. A minimum separation of the bodies is specified for inclusion in this neighbour list. Information on the speed of the particles and the minimum separation distance can be used to estimate the minimum duration of validity of the list. This means that it is impossible for any collision to occur between entities that are not in the neighbour list. In this way only these pairs need to be checked at each timestep creating considerable savings. When this period of validity is exceeded, the process cannot guarantee that entities will not overlap without a contact first being detected. The neighbour list then needs to be updated. Typically this is done by performing a new search.
A DEM process may therefore use two nested loops. The outer loop controls the timing of the search and the construction/or update of the neighbour list. For each incarnation of the neighbour list, the inner loop causes a series of timesteps to be performed, marching the system forward in time using the list to supply pairs to be checked for contact and for any subsequent calculation for those that are. The forces are then summed and the particles and object motions are integrated.
The principles of such a DEM simulation process are illustrated for clarity in FIG. 1 for a given simulation environment, having a mixture of particles with predetermined attributes as listed above, and having predetermined boundary objects (such as the walls of a container, for example).
First, the initial positions (including cell locations), orientations, velocities and spin for each of the particles in the simulation are determined for time t=0 at step 102. At step 104, the neighbour list is determined as discussed above to include pairs of entities (particles or boundary objects). At step 106, pairs of entities in the neighbour list are considered to determine whether they are in contact with each at time t. For each pair of entities in contact with one another or, the net forces on the entities, and their torques are calculated at step 108 and integrated at step 110 to determine new positions, orientations, velocities and spins for the entities concerned.
At step 112, time, t, is incremented and, at step 114, t is compared with the period of validity of the neighbour list to determine whether that period has expired. If it has, the process returns to step 104 and the search for potentially interacting entities is repeated; if it has not, the process returns to step 106 and steps 106 to 114 are repeated.
For the sake of completeness it is noted that, in addition to the collision forces calculated for the DEM method, additional forces can be included to represent the influence of other factors, such as interstitial slurries, inter-particle cohesion, electromagnetic effects and other body forces. These forces are added to the collision forces to give net forces on the particles which are used in the integration of Newton's equation to predict the motion of the particles and objects.
Furthermore, in addition, to determining the motion of the particles and objects, the DEM method also involves generating outputs based on the motion which are used for understanding and describing the processes being modelled. These include, but are not limited to prediction of energy flows and dissipation, wear on the boundary objects, breakage and growth of the particles and information on mixing and segregation.
In the DEM method, the particles can change in number and their attributes can change with time. These can result from flow into and out of the computational domain, breakage of particles, processes to destroy and create particles, agglomeration, erosion and shape change.
Many search methods, performed at step 104 of FIG. 1, have been developed and most rely on the use of a grid to perform localisation of the particles and parts of boundary objects. Some rely on ordered lists, some on tree sub-division and some others on building meshes (such as with Voronoi tessellation). Typically, non-grid methods are not of order N, but are often of order N log(N) making them unsuited to large scale DEM simulation.
Grid based searches use a regular grid of square (in two dimensions) and cubic cells (in three dimensions) which are overlayed over the particles and objects. These entities (particles and object components) are mapped into the cells of the grid and stored in some form of structure for each grid cell. Neighbours of particles in a given cell can be found by checking for entities that lie within neighbouring cells. Since the neighbour cells of a known cell are simple to identify, the grid localises the search from the entire solution space to a known and identified subset whose contents are then easily identified.
The main decision to make when using a grid is the size of the cells to be used and the number of cells around any cell that need to be checked for neighbours. There are three main options:    (1) The cell size is small enough so that the smallest particle in the simulation only just fits into the cell. This guarantees that there can only be one particle per cell and so storage of the particles in a cell is very simple. When finding neighbours of a particle, the search has to be performed by checking enough cells on all sides so that any cell containing the largest particle in the system which is close enough to be a neighbour can be detected. The physical distance that needs to be searched is controlled by the sum of the entity sizes (ie. sum of particle radii for spherical particles). The weakness of this approach is that as the particle size ratio or degree of poly-dispersity (meaning the ratio between the largest and the smallest particles) increases, the number of cells that cover this physical search space becomes very large. For example with a 100:1 size ratio, the largest particle will cover around 100 cells in each direction. The worst case is two proximate large particles. This requires searching a cube of cells of size 3003 (or 27 million cells). This has to be repeated for finding the neighbours of each cell in the grid. The cost of all these checks becomes larger with increasing size ratio and search performance degrades substantially.    (2) The cell size is just large enough to cover the largest particle. In this case, the physical search space is always covered by a cube of 3 cells in each direction—one cell on each side of the base cell whose neighbours are being searched in each of the three dimensions (ie 27 cells). So the identification of neighbour cells is very simple and fast. This method though means that there can be many particles within each cell. Some form of storage structure such as a list, link list or hash table is required to store the identity of particles and object components within each cell. Looking these up adds to the expense of this search variant. This approach works well for modest poly-dispersity but degrades strongly as the particle size ratio increases. The principle reason for this is that the grid localises the search to specific cells very effectively, but the checking of pairs within each cell and between cells is pair-wise for all combinations of entities within these cells. This process is of order M1×M2, where M1 and M2 are the number of entities within the two cells being considered. As the particle size range increases, the number of particles within each cell increases (cubically). Since fine particles tend to segregate in most granular flows, the worst case is a cell full of the finest (smallest) particles being checked against a similar cell. For example, for a size ratio of 100:1, there could be 106 fine particles per cell, so the cost of checks between adjacent cells is of the order of 1012 checks per pair of cells or 1013 per central cell. This cost rapidly becomes extremely high and algorithm performance is seriously degraded.    (3) A cell size in between that of options 1 and 2. This has the weaknesses of both methods in that it needs to store multiple particles per cell and search out a size distribution dependent number of cells around each central cell. The algorithm can degrade both because of chocking large numbers of empty cells or cells with the same large particle or because of concentrations of fine particles leading to large numbers of pair checks between cells.
All variants of the grid search degrade sharply with increasing particle size range. Typically size ratios of 10:1 start to become more expensive than for relatively mono-sized particles. Size ratios of 100:1 tend to become too expensive to simulate and so simulations are restricted to moderate size ranges.
It is desired to address the above or at least provide a useful alternative.